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Abstract 

We present a method for reconstructing two-dimensional velocity fields 
at specified length scales using observational data from tracer particles in 
a flow, without the need for interpolation or smoothing. The algorithm, 
adapted from techniques proposed for oceanography, involves a least- 
squares projection of the measurements onto a set of two-dimensional, 
incompressible basis modes with known length scales. Those modes are 
constructed from components of the velocity potential function, which 
accounts for inflow and outflow at the open boundaries of the measure- 
ment region; and components of the streamfunction, which accounts for 
the remainder of the flow. All calculations are evaluated at particle lo- 
cations, without interpolation onto an arbitrary grid. Since the modes 
have a well-defined length scales, scale-local flow properties are available 
directly. The technique eliminates outlier particles automatically and re- 
duces the apparent compressibility of the data. Moreover the technique 
can be used to produce spatial power spectra and to evaluate the spatial 
effects of open boundaries; it also holds promise for direct calculation of 
scale-to-scale transfer of enstrophy and energy. 



1 Introduction 

Over the past two decades, particle-based video measurements have grown ro- 
bust enough to become standard tools in experimental fluid mechanics. The 
powerful techniques of Particle Image Velocimetry (PIV) and its cousin Particle 
Tracking Velocimetry (PTV) are widely used to obtain Eulerian flow informa- 
tion Tropea et al. 2007], and under proper conditions can be very accurate. 
More general particle-tracking methods can be used to obtain Lagrangian infor- 
mation, allowing the measurement of, for example, the statistics of the velocity 



gradient along particle trajectories Liithi et al. 2005 , the so-called Lagrangian 



Coherent Structures that contain information about the stretching of material 
hues Voth et al. 2002 Mathur et al. 2007| , and the motion of topological 



singularities in the flow Ouellette and Gollub 2007 



All of these techniques involve approximating the flow velocity field — a vec- 
tor function of continuous arguments — from measurements at a discrete set of 
locations. The full field is then usually obtained by interpolation between the 
points, generally either linearly or using cubic splines. Such an approach is of- 
ten reasonable, and is widely applied in computational fluid dynamics (CFD) 
models, where the equations are solved on a discrete mesh. Unlike in CFD, 
however, these experimental velocimetry methods have no control over the lo- 
cations of the mesh points: velocities are known only at the locations of the 
particles, which are randomly positioned throughout the flow. In particular, 
experimentalists cannot employ techniques such as adaptive mesh refinement 
in regions where the flow changes rapidly in space. Simple interpolation, then, 
can introduce significant error into the resulting flow fields. To mitigate such 
effects, experimentalists often use some kind of smoothing filter on the data; 
the choice of filter type and parameters, however, can be more an art than a 
science. We seek instead a data-processing scheme that imposes only a well de- 
fined set of physical constraints on the measurements with as little averaging or 
smoothing as possible, so that the data itself determines the spatial resolution 
and coherence. 

Here, we present an alternate methodology for reconstructing the velocity 
field that has a sound mathematical foundation and requires no ad hoc smooth- 



ing. Our technique borrows tools employed in oceanography Lynch 1989 Chu 



et al. 2003 Lekien et al. 2004 , where researchers have long faced the challenge 



of systematically constructing velocity field estimates from poorly sampled data. 
As we describe in detail below, we first choose a set of basis functions that are 
appropriate for the measured data, imposing physical constraints such as in- 
compressibility on the form of these functions to control the properties of the 
resulting reconstructed velocity field. The data are then projected onto the ba- 
sis functions in a least-square sense Press et al. 2007 ; no additional smoothing 
is required. 

By combining the powerful techniques developed by oceanographers for han- 
dling sparse data sets with the high spatial resolution provided by modern op- 
tical fiow diagnostics, our reconstruction method gives us access to remarkably 
clean and well controlled data sets. It also gives direct access to the scale-local 
field properties: since each basis function has a well-defined length scale, we can 
easily project onto a subset of those modes to construct spatially resolved fields 
with known spectral content (low-pass, high-pass, or band-pass). This type of 
information has recently been exploited to gain a detailed understanding of the 
spatial structures responsible for scale-to-scale enstrophy and energy transfer 

but is typically obtained using more 



Rivera et al. 2003 Bruneau et al. 2007 



complex and less direct methods such as filter-space techniques Rivera et al. 
or wavelet filters 



2003 



Bruneau et al. 2007 



We have developed and tested our reconstruction method using data from 
a quasi-two-dimensional, electromagnetically forced, shallow, stratified flow. In 



2 



^ we describe the flow in detail and discuss the particle-tracking system we 
use to make measurements. In J|3j we introduce the reconstruction technique, 
including the theoretical background for the method and some details of its ac- 
tual implementation. Of particular interest is the choice of basis modes, and we 
discuss the relative strengths and weaknesses of both Fourier modes and eigen- 
modes calculated numerically. In ^ we give examples of how our reconstruction 
scheme can be used to obtain spatially resolved scale-local measurements. Fi- 
nally, in SjSjwe summarize our results and give some directions for future work. 



2 Experimental Methods 
2.1 Flow 

We have designed and constructed an experimental apparatus, sketched in fig- 
ure [T] intended to serve as a physical model of two-dimensional fluid flows. 
Essentially a shallow pan of stably-stratifled salt water, it uses electromagnetic 
forcing and resembles previous devices Tabeling et al. 19911 'Rothstein et al.l 



1999, Solomon and Mezic,,2003tpercx et al.i,2003, Rivera and Eckc, 2005,, Rossi 



et al. 2006 , though it is larger than many. Its test bed has a smooth, glass floor 
86 cm x 86 cm x 3.2 mm thick. We coat the floor with a hydrophobic wax to 
reduce friction, cover it with deionized water, and then inject heavier salt water 
from below such that the stratiflcation is gravitationally stable (denser salt water 
below lighter freshwater). During the experiments discussed below, the saltwa- 
ter solution was NaCl in water, 16% by mass, with a density of p = 1116 'kg/vo? 
(speciflc gravity 1.116) and a kinematic viscosity of = 1.24 x 10~^ m^/s. 

Below the floor of the test bed lies a square array of 34 x 34 neodymium-iron- 
boron (NdFeB) grade N52 magnets, spaced 2.54 cm on center, with alternating 
polarity. The magnets are cylindrical, with a diameter of 12.7 mm and a thick- 
ness of 3.2 mm. Each magnet has a fleld that decays exponentially with the 
axial distance from its surface (see figure[2|, with a peak value near 0.3 T and a 
decay length of C = 4.70 mm. A pair of bar electrodes are mounted on opposite 
ends of the test bed and connected to a BK Precision 1670A power supply, al- 
lowing an electric current to be passed through the saltwater layer. That current 
interacts with the magnetic flelds via the Lorentz force F — J x B/p, where 
J is the current density and B is the magnetic fleld, to produce bulk motion 
in the salt water. For all experiments discussed below, the forcing current was 
steady in time. At low current density J — |J|, the forcing is gentle and the 
flow generated is a steady lattice of vortices of (spatially) alternating direction. 
At higher currents, stronger forcing leads to time-dependent flows, loss of lattice 
symmetry, spatiotemporal chaos, and perhaps two-dimensional turbulence. We 
characterize the flow by its Reynolds number 

^^^w™^^ (1) 



where Lf — 2.54 cm is the forcing scale and u^ms ~ \J {u^ + v^) is the root- 
mean-square velocity, with (•) signifying an average over all particles in each 
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Figtire 1: The experimental apparatus, as seen (a) empty from above, and (b) 
filled in cross-section. Both sketches are partial; the full test bed is square, 
with walls on all sides. Electrical current flows from left to right through the 
saltwater layer. 
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Figure 2: Axial variation of the field of an example magnet, as measured with 
a Gaussmeter. The plotted exponential fit has Bq — 0.28 T and = 4.70 mm. 



frame. Here and throughout we represent the two-dimensional flow as m = 
ux + vy, adopting Cartesian coordinates {x,y,z) with unit vectors {x,y,z), 
respectively, x being the direction of the forcing current and z being the out- 
of-plane direction. Flows with 50 < Re < 250 are accessible with our current 
experimental setup. 

To the saltwater layer we add tracer particles whose motion indicates that 
of the fluid. The particles are fluorescent polystyrene spheres, 51 fiui in diam- 
eter. Their specific gravity is 1.05, intermediate between that of the saltwater 
layer and the freshwater layer. Thus they are trapped at the interface be- 
tween the layers. Due to the miscibility of the two fluid layers, there are no 
long-range surface-tension forces between the particles. In contrast, if the parti- 
cles are allowed to float to an interface between salt water and air, they quickly 



clump Vella and Mahadevan 20051 and become poor tracers, with a much larger 



effective diameter. The particles absorb most strongly in the blue (468 nm) and 
emit most strongly in the green (508 nm). We illuminate them with blue light 
emitting diodes (LEDs) whose luminosity peaks at 470 nm. We subsequently 
image the particles using an IDT MotionPro M5 camera, first filtering optically 
to attenuate wavelengths below 520 nm, thereby reducing blue glare much more 
than the green light emitted by the particles. The camera has a CMOS sensor 
with pixel dimensions 2320 x 1728 and can record continuously at frame rates 
up to 170 Hz. We trigger each frame externally, using an Agilent 332 10 A func- 
tion generator. The resulting movies are recorded to disk for later processing 
and particle tracking. 

2.2 Particle Tracking and Velocimetry 

To measure the flow dynamics, we identify and follow the particles in recorded 



movies using Lagrangian particle tracking Ouellette et al. 2006 Ouellette 
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2010 . As is done in classical PTV, we locate the positions of individual particles 
and match them in time. Unlike in PTV, however, which typically only uses 
two images to match particles, we keep the full Lagrangian information about 
the long-time particle trajectories. Lagrangian dynamics are a powerful way to 
study fluid flow Toschi and Bodenschatz 2009 ; for the present application, the 
length of our trajectories is useful only in that it allows us to use a numerical 
differentiation scheme that is more robust than simple finite differences. 



Our tracking algorithm has been described in detail elsewhere [Ouellette 
et al. 2006 [Ouellette 2010 , and so here we only briefly summarize the major 



steps. Particles are located by using local maxima in image intensity (above 
a small threshold) and the pixels immediately adjacent to the maxima. We 
obtain the particle centers with a resolution of roughly 0.1 pixels (13.6 /zm in 
the experiments described herein) by fitting one-dimensional Gaussians to the 
horizontal and vertical intensity profiles. Once the particle locations have been 
found, they are tracked using a predictive three-frame best-estimate algorithm. 
For each partially constructed trajectory, the expected position of the particle 
at the next time step is estimated using simple kinematics. The measured 
particle position that comes closest to this estimate is then chosen to continue 
the trajectory. Small temporal gaps (due to particle drop-out) are bridged via 
extrapolation. The entire process has been shown to be robust under a wide 
range of flow conditions, and has been parallelized so as to be computationally 
efficient. For the data reported here, we tracked Np ~ 20 000 particles per 
frame. We subsequently differentiate the trajectories temporally by convolving 
the measured trajectories with a Gaussian smoothing and differentiating kernel 



Mordant et al. 2004 Ouellette et al. 2007 , yielding a time series of positions 



and velocities for each tracked particle. 



3 Velocity Reconstruction 

After performing Lagrangian particle tracking and differentiating the resulting 
trajectories, we are left with a collection of particle tracks that, taken together, 
approximate the velocity field of the fluid. The state of the art in imaging tech- 
nology and computer processing is such that we can produce data sets whose 
size and level of detail approaches that of numerical simulation. Other charac- 
teristics of our data sets, however, are unlike numerical data — in many ways 
they resemble the observational data of oceanographers. First, our velocity 
measurements do not lie on a regular grid, but lie at locations that change from 
frame to frame as the tracer particles follow the flow. The spatial distribution of 
the tracers also tends to be inhomogeneous. Second, outlier points unavoidably 
appear in the data. A small fraction are due to physical imperfections in the 
experimental apparatus, such as dust particles floating on the water surface. 
More arise through the course of tracking particles and differentiating tracks. 
Whatever their source, clearly non-physical outlier points require correction. 
Third, a small but finite compressibility appears in our velocity measurements. 
Since the typical flow velocities are on the order of 1 cm/s and the speed of 
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sound in water is 1484 m/s, we do not expect physical compressibility; rather, 
the apparent compressibility is due partly to vertical flow and largely to nu- 
merical noise inherent in spatial gradients calculated from irregular samples at 
finite resolution. These three characteristics, common to any experimental data 
gathered by PTV or by Lagrangian particle tracking, present challenges and 
highlight the need for careful post-processing. 

Perhaps there are other ways to sidestep these three challenges. PIV mea- 
sures fluid velocities on a regular grid, avoiding the first challenge, but PIV 
is necessarily Eulerian. Interpolating the velocity field onto a regular grid is 
straightforward but requires smoothing with parameters that are typically em- 
pirical and semi-arbitrary. Perhaps worse, interpolation can yield more appar- 
ent information than has actually been measured. Adjusting the algorithmic 
parameters used in identifying and tracking particles can reduce the number 
of outliers substantially, addressing the second challenge. Outliers cannot be 
eliminated entirely, however, and reducing their count means an even larger 
reduction in the count of "good" particles being tracked — valuable data is 



lost Ouellette et al. 2006 Ouellette 2010 . A variety of experimental tech 



niques have been developed to reduce vertical flow in devices like ours, and we 
employ many such techniques to address the third challenge. Interpolating and 
smoothing might lead to spatial gradients with less numerical noise and thus a 
lower apparent compressibility, but the choice of smoothing and interpolation 
parameters remains empirical at best. 

The post-processing technique presented below addresses these challenges 
while producing a final data set that more directly resembles the experimental 
measurements, relying on few physical or numerical assumptions. Moreover, 
it allows for multi-scale reconstruction of the velocity field at specified length 
scales. The remainder of this section details the technique, first outlining the 
construction of basis modes for the velocity components from a velocity po- 
tential function and a streamfunction, then describing numerical techniques for 
calculating spatial gradients without interpolation. Continuing, we present the 
projection algorithm, show how to account for open boundaries with velocity 
potential modes, address the choice of basis for the streamfunction modes, and 
close with a concise summary of the algorithm. 



3.1 Velocity Basis 



Two-dimensional, incompressible vector fields can always be represented in 
terms of a streamfunction, a mathematical fact that finds frequent use in the 
study of fluids. Its validity arises from the Helmholtz-Hodge theorem: any vec- 
tor field that vanishes at its boundary (possibly at infinity) can be uniquely de 
composed into an incompressible component and an irrotational component Ar^ 
fken and Weber 2001 . In two dimensions. 



V$ - i X V*, 



(2) 
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where 



dx 



y 



dy' 



(3) 



In the specific case where m is a velocity field, $ is known as the velocity 
potential and ^' as the streamfunction. The first term on the right-hand side of 
equation [2] is irrotational in two dimensions, and the second is solenoidal in two 
dimensions, that is. 



V X (V$) = 0, 
V-(ixVI') = 0. 

Alternatively, the streamfunction may be understood in terms of toroidal 
and poloidal components. Any three-dimensional, incompressible vector field 
can be decomposed uniquely into a toroidal and a poloidal component, as is 



commonly done in geophysics, for example by BuUard and Gellman 1954 . The 



symmetry axis of the toroidal and poloidal components may be chosen a priori 
and should correspond to the symmetry of the system being considered; with 
in-plane flow in mind, we can choose the out-of-plane vector z and write the 



decomposition as Chu et al. 2003 



M = V3 X ivj/ + V3 X V3 X zx, 



(4) 



where \l/ and x are the toroidal and poloidal scalar functions, respectively, and 



V3 = * 



d_ 

dx 



y 



d_ 
dy 



dz 



(5) 



In the special case of incompressible, two-dimensional flow, x = and ^' again 
becomes the streamfunction. 

These mathematical identities immediately allow for sophisticated manip- 
ulation of velocity fields that are already incompressible and two-dimensional. 
They also allow for conditioning of imperfect, real- world data by representing 
those measurements in terms of a two-dimensional, incompressible basis — be- 
cause such a basis can be built by decomposing the velocity potential function 
and streamfunction. Following Chu et al. 2003 and Lekien et al. 2004 , we 



write the velocity potential function as 



^{x,y) ^^aj(j)j{x,y), 



(6) 



with some coefficients aj, for some number A^f, of orthogonal basis functions 
4>j{x,y), which we will call the "boundary modes" and usually abbreviate cjjj. 
Likewise we write the streamfunction as 



Ni 

E 

k=l 



I3ki^k{x,y), 



(7) 
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with some coefficients for some number Ni of orthogonal basis functions 
tjjkix,y), which we will call the "interior modes" and usually abbreviate ipk- 
Choosing Nt, ~ Ni = oo allows exact representation of any velocity potential 
function and any streamfunction, but in numerical calculations it is necessary 
to truncate the series by choosing Nh and Ni finite. The choice of Nf, and Ni, 
as well as the exact forms of the boundary modes and the interior modes ipk-, 
will be specified below. 

As mentioned above, one hypothesis of the Helmholtz-Hodge theorem (equa- 
tion [2]) requires that the amplitude of the vector field be zero at its boundary. 
However, we typically observe particle velocities in a central sub-region of our 
experimental apparatus, such that the rigid sidewalls are not in view, and flow 
through the boundaries of the measurement region is not only possible but 
ever-present. Thus the boundary requirement is not satisfied. This situation 
is common in experimental and observational studies, where measurements are 
often recorded in a region smaller than the fluid container or oceanic coastline. 
Mathematically, such boundaries are said to be "open." With open bound- 
aries, decomposition into incompressible and irrotational components remains 



possible, but is no longer unique Lynch 1989 . Rather, the results of de 



composition depend on the boundary conditions, which can be chosen freely. 



Following Lekien et al. 2004 , we choose the boundary condition 

*lr = 0, (8) 

implying also i/'feir ~ "where F represents the boundary. In words, all inflow 
and outflow at the boundary is accounted by the velocity potential function $, 
not the streamfunction ^. 

With equations [2] and [8) the Helmholtz-Hodge decomposition of the velocity 
is specified uniquely, and with equations [7] and [6j the streamfunction and velocity 
potential function have been constructed. However, in experiments we do not 
measure the streamfunction or velocity potential function, so we cannot project 
our data onto them directly. Rather we build basis sets for the two horizontal 
components of the velocity by using equation [2] with each basis mode, yielding 

i=i k=i ^ 



j=i fc=i 



Hence we have constructed incompressible, two-dimensional basis sets onto 
which measurements can be projected. Since the velocity components u and 
V are constructed from the same velocity potential function $ and same stream- 



function ^I*, the coefficients aj and Pk are the same for both; see equation 12 
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3.2 Gradients at Particle Locations 

As equation [9] makes explicit, constructing the velocity basis from the velocity 
potential function basis 4>j and the streamfunction basis ipk requires calculating 
their spatial gradients. If (jjj and V'fc were known analytically, their gradients 
could be calculated exactly. On the other hand, if and tpk are constructed 
numerically, their gradients must be approximated. As emphasized above, our 
intention is to condition experimental measurements and gain access to their 
properties at various length scales, while maintaining as much fidelity to the 
original measurements as possible. We would like to minimize the number of 
post-processing parameters introduced, avoiding interpolation and smoothing. 
Our method for numerical calculation of gradients must follow suit. 

Our observations, however, are not made at spatial locations that are reg- 
ular or predictable — we measure the fluid velocity at the particle locations 
(which are effectively random) and nowhere else, as figure [S] shows. Thus grid- 
based methods for calculating gradients, such as simple finite differences, would 
require interpolation onto a regular grid and are inappropriate here. When it 
becomes necessary to calculate numerical gradients for use in equation |9] or in 
equations below, we do so at the particle locations via a finite-element technique. 
Typically, finite-element methods are used to solve partial differential equations 
on the vertices of a triangular "mesh" that is constructed a priori and often 
refined iteratively to produce sufficiently low numerical error. Our approach 
differs: the mesh is not constructed nor refined a priori, but is specified by the 
particles themselves, whose locations are the vertices of the mesh. Its edges are 
constructed via a Delaunay triangulation, after which we remove edge triangles 
whose aspect ratios are too small (below 0.1) and who are therefore prone to 
numerical instability. Consequently, the velocity bases given by equation [9j like 
all other quantities, are known at the particle locations, and neither interpola- 
tion nor smoothing is necessary. A new mesh like the one shown in figure |3] is 
constructed for each frame of data, since particles move between frames. 



3.3 Least-Squares Projection 

Once incompressible, two-dimensional bases have been constructed at the parti- 
cle locations for both components of the velocity, the next task is to determine 
the coefficients aj and /3k in equation [9j They can be calculated by projecting 
measurements from tracked particles onto the bases using a linear least-squares 
method. For notational simplicity, we define the full set of basis functions (in- 
cluding contributions from both the velocity potential function and the stream- 
function) as 

n - i5Ju^*|, (10) 



Vl = 





H 






H 


I dy J 



dy 

dipk 
dx 
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Figure 3: A portion of an example particle mesh as seen from above. Particle 
locations are indicated by black dots much larger than the actual size of the 
particles. The mesh itself, computed from the particle locations via a Delaunay 
triangulation, is shown in grey. 



and the full set of coefficients as 7/ — {ctj} U {/3/c}- Then equation [9] can be 
rewritten as 



1=1 
1=1 

The total squared error of the projection is then the square of the difference 
between equation evaluated at the particle locations, and the measurements 
themselves. Minimizing that error is mathematically identical to maximizing 
the probability that the coefficients 7; (or equivalently, aj and /3k) correctly 



model the measured data, assuming Gaussian errors. It can be shown Press 



et al. 2007 that the error is minimized by the coefficients 7; that satisfy the 



matrix equations 

E E K^,p+«;pt',p)7; = EK"r"'+^^p<'"0' (12) 
p=i 1=1 p=i 

known as the "normal equations" of the linear, least-squares fit. Here p indexes 
the Np particles present in the frame, and all sums are written explicitly. The 



notation uip signifies the value of the basis mode ui, as defined in equation 10 
evaluated at the location of particle p; Ujp, vip, and Vjp have corresponding 
meanings. Likewise and w™^'^'* are the components of the velocity as 

measured at particle p. 



To solve equation 12 we use singular value decomposition, a mathematical 



operation that avoids numerical issues that arise due to nearly singular matrices 
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and due to basis modes that are nearly degenerate at the particular locations 



where measurements have been recorded Press et al. 2007 . Once the coeffi- 
cients 7/ are known, the velocity components u and v can be reconstructed by 
simply summing over the basis modes, that is, by using equation 1 11 1 directly. 

3.4 Boundary Modes 

So far we have required only that the basis functions 4'j and ipk be orthogonal, 
without specifying them exactly; we now return to the question of choosing 
those bases, first discussing the boundary modes then the interior modes 

To construct a basis for the velocity potential function, we follow the method 
of 'Lekien et al. 2004 . In an incompressible, two-dimensional flow with closed 



boundaries, $ = identically. We retain the velocity potential function, how- 
ever, to account for inflow and outflow through the open measurement bound- 
aries — a choice that led to equation [Sj Thus the boundary modes are deter- 
mined completely by the boundary conditions. We denote the local outflow at 
the boundary F with the function G'(t) = n ■ u, where n is the unit vector 
locally normal to the boundary and pointing outward, and the coordinate r 
gives the distance around the boundary, measured from some arbitrary starting 
point. The local outflow can be decomposed in Fourier modes as 

G{r)=Y,^^,gJ'^y (13) 

j=i \ / 

with some coefficients fij, for some number Nb of modes, where Li, is the total 
length of the boundary (so < t < Lb) and 

, ^ J sinjx, j even, , , 

9j[^) - I ^ j^^^;, j odd. ^ ' 

We choose the basis this way so that each gj is continuous over the boundary 
and ^ ^ 

G(r)dT = ^ ' g, (^^^ dr - 0, (15) 

as required by incompressibility. 

We construct the basis modes for the velocity potential function using so- 
lutions of Laplace's equation, with the outflow modes gj{'KT /Lb) providing the 
Neumann boundary condition: 

V2(/)j = (16) 



^ , , , 7rr 

n ■ V 



That the modes (pj must satisfy Laplace's equation is also a requirement of in- 
compressibility, made clear by taking the divergence of both sides of equation [2j 
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Figure 4: Example basis modes, (a) Three boundary modes (jjj. A dot is drawn 
at the location of each tracked particle, much larger than the actual particle, 
with its shade representing the value of the mode at that location. Dark shades 
signify positive values and light shades, negative. The amplitude is arbitrary, 
and the white background shows through in regions where few particles are 
present. The particles are seen from above, and the forcing current flows from 
left to right. The value of j is indicated in each case, (b) Three interior modes 
tpk, drawn in the same way. The values of (n^, Uy) are indicated in each case; see 



equation 19 The velocity bases {ui,vi) are constructed from spatial gradients 



of (pj and ipk according to equation nO] 



Three example basis boundary modes, constructed this way, are shown in fig- 
ure |4] To set TV;,, the total number of boundary modes, we choose a minimum 
length scale Lmm and retain all modes whose scales Lj, satisfy 



3.5 Interior Modes 

To construct a basis for the streamfunction, we use the eigenf unctions of the 
Laplacian, satisfying 

V^ijk = AfcVfc, (18) 

where is the scalar eigenvalue of the eigenfunction ipj.. The eigenfunctions 
of the Laplacian are used in many situations where an orthogonal, complete 
set is required; their exact form varies with the geometry and conditions at 
the boundary. Some examples include spherical harmonics, vector spherical 
harmonics (toroidal and poloidal components, as in equation ffl), and Bessel 
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functions. For a scalar function (such as the streamfunction ^{x,y)) in a two- 
dimensional, rectangular region, the eigenfunctions of the Laplacian are two- 
dimensional Fourier modes (products of sines and cosines). 

For geometries in which eigenfunctions of the Laplacian are not known ana- 
lytically, they may be calculated numerically. Such an approach allows analysis 
of data sets with boundaries of arbitrary geometry, making it valuable for the ir- 
regular observation regions common in oceanography. Its drawback, however, is 
substantial computational expense. In reconstructing observed velocity fields, it 
is necessary to choose between analytical and numerical computation of eigen- 
functions (or "eigenmodes" ) , weighing the strengths and weaknesses of each. 
Below we compare and contrast reconstructions of the same data set using both 
techniques. 

The two-dimensional Fourier basis for the streamfunction takes the form 

. / n^Trx\ . ( nyTTy\ 
^, = sm(^— jsm(^ — j, (19) 

where and Ly are the sizes of the rectangular measurement region in the x 
and y directions, respectively, and k stands for the combination of and rij,, 
both of which are positive integers. Three such modes are shown in figure J4| 



Note that the eigenvalue equation ( 18 1 and the boundary condition (equation^ 
are both satisfied automatically; cosine modes, on the other hand, would not 
satisfy the boundary condition. Fourier modes have well-defined length scales in 
both dimensions, growing smaller as the integers Ux and Uy increase. To set iV^, 
the total number of interior modes, we employ a criterion similar to|17[ retaining 



all modes whose scales are larger than the minimum in both dimensions 

2L 

n 



''X 



2Ly 



^ Lmim (20) 

Uy 

With the criteria [20] and [17) we have all the information necessary to use 



the boundary basis obtained from equation 16 and the interior basis defined 
in equation [19] to reconstruct a velocity field. The effects of varying Lmin are 
summarized in figure [5] Smaller values of Lmm increase Ni and iVf,, captur- 
ing more of the information recorded in the original measurements (which we 
roughly estimate using kinetic energy), but also requiring longer computation 
times and additional storage. One noteworthy feature is the large jump in ki- 
netic energy content near the forcing scale. Consistent with intuition, this jump 
suggests that in order to capture important features in the fiow, one should 
choose Lmin < Lf. Calculation time is plotted on a semi-logarithmic scale — 
its dependence on is steeper than exponential. The calculation time also 

depends strongly on the particle count Np (not plotted). 

Alternatively, a basis for the interior modes can be constructed by solving 
equat ion [T8| numerically. We do so using a finite-element technique similar to the 
one described above in order to obtain numerical basis functions at the particle 



14 




Figure 5: Characteristics of example bases built with varying values of the 
minimum feature size L„yim including (a) the total number of basis modes (iVj 
interior modes and boundary modes), (b) the fraction of kinetic energy 
captured by the reconstruction, and (c) the time required to build the bases 
and calculate the reconstruction for one frame, composed of 19 000 particles. 
Here {KE) is the mean particle kinetic energy after reconstruction, {KEq) is the 
mean particle kinetic energy as measured (before reconstruction), and Re = 100. 
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locations, without interpolation or smoothing. As in the case of Fourier modes, 
it is necessary to decide how many modes Ni are required, which can be done 
by estimating the length scale of each numerical eigenmode using dimensional 



arguments Lekien et al. 2004 . We set N.i by retaining all modes with 



Ll^l^ > Ljnin, (21) 

where Ai is the lowest eigenvalue and Li is the boundary length scale (we use 



Li = L^). With the criteria 17 and 21 we can use the boundary basis obtained 



from equation |16| and the interior basis obtained from equation[l8]to reconstruct 
the velocity field. The choice of Lmin has similar effects as with the Fourier 
basis described above, as shown in figure [Sj The kinetic energy fraction again 
increases sharply near the forcing scale, as one might expect. 

We are now in a position to compare and contrast Fourier modes and nu- 
merical eigenmodes as the choice for the interior modes. As stated above, in 
a rectangular measurement region, the Fourier modes satisfy equation [18] ex- 
actly, and we would expect no difference between the two bases. But even 
though our camera has a rectangular field of view, our measurement region is 
not rectangular — rather, its edges are set by the locations of the outermost 
particles, as shown in figure [3j Thus the Fourier basis and the eigenmode basis 
differ. Figure [5] shows that the value of Lmin where kinetic energy suddenly 
increases is slightly higher for eigenmodes than for Fourier modes, suggesting 
that eigenmodes reproduce the flow more efficiently in terms of £mm- In terms 
of calculation time, however, they are much less efficient. Even though more 
Fourier modes are required to reconstruct a velocity field with sufficient energy 
content, they require much less computation time (by a factor of five in this flow 
with 20 000 particles, if we wish to capture 30% of the measured kinetic energy). 
Thus for the nearly-rectangular boundaries common in our experiments, Fourier 
modes form a more attractive basis than eigenmodes. In measurement regions 
whose boundaries are more irregular, the fidelity of the Fourier basis decreases, 
making the eigenmode basis more attractive. In a circular measurement region, 
Bessel modes would be the most attractive basis set, though we do not consider 
them here. 

Another comparison between the Fourier basis and the eigenmode basis is 
provided in figure |6) which shows the vorticity u; = V x tt of a single frame as 
measured and as reconstructed with each basis. By eye, both reconstructions 
are similar, with noticeable noise reduction when compared to the measured 
data, though Fourier modes are computationally less expensive. Henceforth we 
shall use Fourier modes with Lmin = Lf/A. 

The kinetic energy fraction plotted in figure [s] does not reach 100% and 
perhaps deserves comment. Where is the kinetic energy going, and would inter- 
polation preserve more of it? Particle tracking, like any measurement technique, 
involves systematic uncertainty which in this case tends to increase the kinetic 



energy artificially. Ouellette et al. 2007 have estimated the uncertainty asso 
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Figure 6: An example vorticity field seen from above (a) as measured, (b) after 
reconstruction using boundary modes and interior Fourier modes, and (c) after 
reconstruction using boundary modes and interior cigenmodes. A dot is drawn 
at the location of each tracked particle, much larger than the actual particle, 
with its shade representing the value of the vorticity at that location. The white 
background shows through in regions where few particles are present. For both 
projections, Lj^in = 0.27i/. Here Re = 100, and the forcing current flows from 
left to right. 
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ciated with tracking and differentiating the tracks as 



= (22) 

where ^^2: (2^)^^ is the velocity variance due to uncertainty, 2T is the integra- 
tion time used when convolvingthe measured trajectories with the smoothing 
and differentiating kernel (see ^2|, and (e^) is the variance of the uncertainty 
in locating particles spatially. In our measurements, T = 1.5 frames and {ex) 

is about 0.1 pixels, yielding {^^3. (T)^^ — 0.0063 pixels^ /frame^, which corre- 
sponds to a standard deviation about 3.3% of the root-mean-square velocity at 
Re = 100. Additional uncertainty arises because of the presence of physical 
imperfections in the experimental apparatus, such as dust or tracer particles 
floating on the water surface. Ideally, both types of uncertainty would be elim- 
inated in post-processing. Hence we do not expect the kinetic energy fraction 
of any reconstruction to reach 100%. 

One way to examine the effect of reconstruction in more detail is by plot- 
ting probability density functions (PDFs) of velocity. We have reconstructed 
50 frames of measured data {Re — 100) by projection with a Fourier basis for the 
interior modes, using Lmin — Lf We have also reconstructed the same data 
with interpolation and smoothing, using a smoothing length of L//4. Figure [t] 
shows PDFs of the measured data and both reconstructions. The measured 
data shows non-Gaussian tails in its velocity PDF which are removed by either 
reconstruction technique, as one would hope. The kinetic energy fraction is re- 
lated to the variance of the velocity, indicated by the width of each PDF. As the 
figure shows, reconstruction by projection results in a velocity variance much 
closer to the original than does the more common interpolation and smooth- 
ing — thus more of the kinetic energy is preserved. In this example, projection 
preserves about 33% of the kinetic energy, while interpolation and smoothing 
preserves only about 2.3%. Thus the kinetic energy fraction plotted in figure [s] 
would be considerably lower if we used interpolation instead of projection. 

3.6 Summary of the Algorithm 

Readers who are interested in reconstructing other velocity fields using this 
method may find a concise summary useful. We begin by building a parti- 
cle mesh via Delaunay triangulation, identifying the edge particles. Next we 



construct the outflow modes gj{TTT/Lf,) (equation 14|, using equation 17 to de 



termine the mode count. The gj give the boundary conditions for a numerical 



solution of Laplace's equation (equation 16), for which we use finite-element 
techniques, yielding the boundary modes (f)j. We calculate their gradients nu- 
merically to obtain the components of the velocity bases (equation [9]) that arise 
from the boundary modes. For the special case of closed boundaries (without 
inflow or outflow), the boundary modes may be ignored altogether, that is, 
6, = 0. 
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Figure 7: Probability density functions of the velocity component v = u ■ y, 
comparing reconstruction by projection (using the Fourier basis for internal 
modes) to interpolation and smoothing. In this example projection preserves 
33% of the original kinetic energy, while interpolation and smoothing preserves 
2.3%. The measurement encompasses 50 frames, including about 10^ particles, 
at Re = 100. 
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The interior Fourier modes ipk, meanwhile, are given by equation 19 and 



their gradients can be calculated analytically. Thus the components of the ve- 
locity bases (equation [9]) that arise from interior modes can be constructed ex- 
plicitly from trigonometric functions, the mode count being set by equation |20[ 
With the velocity basis now complete, we solve the normal equations of the 
linear, least-squares fit (equation 12 ) numerically using singular value decom- 
position to obtain the coefficients 7;. Finally we reconstruct the velocity field 



by summing its components according to equation 11 



4 Multi-scale Structure 
4.1 Spectral Power 

Reconstructing observed velocity fields using the technique outlined in f|3] can 
be useful to condition a data set, removing outlier particles and reducing the 
apparent compressibility associated with flow in the vertical direction. Because 
it involves projection onto a basis of orthogonal modes with well-defined length 
scales, it also offers another benefit, direct access to the multi-scale structure of 
the flow. For example, instantaneous spatial power spectra can be calculated 
without interpolation or smoothing — the spectral power in each mode is simply 
the square of its coefficient 7;. Producing a one-dimensional spectrum requires 
averaging over one spatial dimension, and a variety of choices of its direction 
are possible, including x, y, and the azimuthal direction. Averaging the spectral 
power in the y direction and sorting the averaged modes according to their x 



wavenumber — 2Tmx/Lx (see equation 19) yields spectra like those shown 
in figure |8] Each is the product of averaging 36 y-direction modes in each of 
997 frames of data, spanning about 1 min in time and including roughly 10^ 
particles. No interpolation, smoothing, or windowing is required. 

Consistent with our expectations, figure [8] shows a peak in power at the 
wavenumber nearest the forcing scale Lf. The spectra are robust in that they 
account for a large number of measurements, but their range of scales is insuffi- 
cient to comment on scaling laws or on the presence or absence of the enstrophy 
and energy cascades predicted by Kraichnan] [1967 . Choosing a smaller value 



for Ljnin would add modes at smaller scales and increase the range, but the 
increase in computational expense is steep (see figure |5| . When a spectrum is 
desired over a wider range of scales, a better method would be to use criteria 
different from equation [20] for choosing the number of basis functions in the x 
and y directions, such as 



TTUx 

Tin,, 



With the dimensionless parameter ^ < 1, this basis would include more scales 
in the x direction without adding more modes in the y direction. Thus the 
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Figure 8: Normalized power speetral density of the veloeity field at various 
Reynolds numbers, calculated using the coefficients /3k of basis modes with well- 
defined length scales. The forcing scale L j is indicated with a vertical line. Here 
(•)^ signifies averaging in the y direction, that is, averaging all modes with the 
same kx- 



21 



computational expense would remain manageable, and additional variation in 
y would have been removed in averaging in any case. Since ^ and Lmin can be 
adjusted, no limits to the spectral range need be introduced in post-processing; 
only the limit imposed by the data itself — the particle spacing — remains. 



4.2 Scale-Local Velocity 

Spectra give some measure of the multi-scale structure of a flow and are common 
in the literature because of the relative ease with which they can be obtained, 
via a variety of methods. Reconstruction by projection onto basis modes gives 
access to much more information, however, including scale-local velocity fields. 



Inspired by techniques used in large-eddy simulation, Rivera et al. 2003 de 



veloped a method for calculating the transfer rates of enstrophy and energy 



directly, using band-pass spatial filters. Later Bruneau et al. 2007 introduced 
a similar technique based on wavelets. We suggest that projecting observational 
data onto basis modes with well-defined length scales can allow for the same 
sort of direct calculation of enstrophy and energy fiux, without the complication 
of filter construction and tuning. 

As a proof of principle, we offer figure [9] It reproduces an example vorticity 
field (the same data as in figure[6| , splitting the fiow into components with length 
scales near the forcing scale, larger than forcing scale, and smaller than the 
forcing scale. Instead of applying band-pass filters, we can select only the part 
of the velocity field relevant to a particular band of length scales by retaining 
only those Fourier modes whose length scales 




L=^.\^ + ^ (23) 



fall within the band. The factor 1/^2 assures that modes with 



Lf (24) 



have length scale L = Lf. Choosing 2L//3 < L < 5Lf as the limits of the 
band near the forcing scale, we find that flow in this range of scales bears a 
strong resemblance to the full flow as shown in flgure|6] One might expect the 
resemblance not only because we have included the forcing scale but also because 
the majority of the energy is captured by this band, as indicated by figure [8j 
Plotting the vorticity at different length scales shows that the majority of the 
enstrophy is also captured by the band near , as the color scales of figure |9] 
make clear. The flow features with scales L < 2L//3 appear to be dominated by 
frequency harmonics of the forcing lattice, though other phenomena are present 
as well. The flow features with scales L > 5Lf are dominated by structures that 
are large in one spatial dimension but small in the other. The power associated 
with circulations that are large in both spatial dimensions is small compared to 
other flow components. 
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Figure 9: Components of the vorticity field shown in figure [6] drawn the same 
way, at length scales (a) near the forcing scale, with 2L//3 < L < 5Lf, (b) 
larger than the forcing scale, with 5Lf < L, and (c) smaller than the forcing 
scale, with L < 2Lf/3. Note that the limits of the color scale are different in 
each image. Here we use a Fourier basis for the interior modes, and the length 
scale L of each mode is defined according to equation [23j 
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4.3 Boundary Effects 

In addition to giving access to scale-local spectral power and velocity (and vor- 
ticity) components, the reconstruction method described in Sjs] allows separation 
of the scales and structures associated with boundary effects from those associ- 
ated with interior flow. For two-dimensional, incompressible flows without out- 
flow at the boundaries, the streamfunction as decomposed in equation [7] with 
the boundary condition [8] comprises a sufficient representation, and we need 
no more components. In the case of a flow with open boundaries, the velocity 
potential modes (equation [6| are included as a necessary and sufRcient account 
of the outflow. Since the flows associated with velocity potential modes 4>j are 
orthogonal to those associated with streamfunction modes V'fe, their motions 
may be examined independently. Moreover the dependence of boundary effects 
on flow characteristics and post-processing parameters can also be quantified. 

For example, the fraction of total kinetic energy associated with the bound- 
ary modes varies with Lmm, as shown in figure [TOj When Lmin is chosen near 
the forcing scale Ly, boundary modes are relatively unimportant, accounting 
for less than 2% of the total kinetic energy. As Lmm increases above the forc- 
ing scale, the boundary modes become more important, likely because interior 
velocity variations at scales near Lf are being replaced with their near-zero 
average values, decreasing the contribution of the interior modes. As L^i^ de- 
creases below the forcing scale, boundary modes again become more important. 
The origin of this effect is less clear, but may relate to under-sampling. In the 
frame of data used to produce figure [TOj the measurement region is bounded by 
183 particles, distributed around a region that is nearly rectangular and has a 
perimeter near 110 cm; thus their mean spacing is 0.6 cm. Though the sampling 
theorem applies only in the case of equal spacing, one would nonetheless expect 
the smallest resolvable scales to be nearly twice that length (1.2 cm), and that 
aliasing effects might distort measurements at smaller scales. Figure [T0| shows 
the relative kinetic energy of the boundary modes beginning to increase at scales 
near 1.2 cm, consistent with this hypothesis. 

The choice of boundaries affects the importance of the boundary modes more 
dramatically than the choice of Lmin- By examining sub- regions of our camera's 
field of view, we observe that boundary modes account for a larger fraction of 
the total kinetic energy when the region is small than when the region is large. 
Since the flow is dominated by a lattice of vortices with well defined length scale, 
boundary modes become unimportant only in the region more than one vortex 
diameter from the boundary. Keeping the vortex scale constant, a decrease in 
the boundary size by some factor C reduces the area where boundary modes 
are unimportant by a factor C^. Thus boundary modes have greater influence 
in smaller measurement regions. 

Changing the location of the measurement region, even while maintaining its 
size, can also affect the structure and scales accounted by the boundary modes. 
Consider the two boundaries sketched in figure [Tlj one drawn as a solid line and 
cutting through vortex centers, the other drawn as a dashed line and passing 
near the separatrices between vortices. Though incompressibility requires that 
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Figure 10: Fraction of the total kinetic energy accounted by the boundary 
modes, as a function of the minimum feature size Lmim m an example frame 
of data with Re — 100. Here {KEi,) is the mean particle kinetic energy in the 
boundary modes, and (KE) is the mean particle kinetic energy after reconstruc- 
tion with all modes. 



the net outflow be zero in both cases, the local inflow and outflow (equation 13 1 
is evidently much larger for the first boundary than for the second. Accordingly, 
the boundary modes are more important. Reconstructing 50 frames of the same 
data with each boundary, we find that in the first case, boundary modes account 
for 18.4 ± 7.4% of the kinetic energy, whereas in the second case, boundary 
modes account for just 8.1 ± 3.0% of the kinetic energy. 

The fraction of total kinetic energy accounted by the boundary modes also 
varies with other parameters. Preliminarily we observe that it increases with 
Reynolds number. Larger Reynolds numbers may drive more of the enstrophy 
to scales smaller than Lmim causing the interior fiow to be damped by averaging 
in the same way we have postulated above, for the case of large L^nin- 

Finally, before leaving the topic of boundary effects, we note that we have 
also implemented a reconstruction algorithm more closely resembling that de- 



scribed by Chu et al. 2003 . Instead of accounting for outflow at open bound- 
aries by using the velocity potential function, it adjusts the streamfunction 
modes directly by solving the eigenvalue problem [18] with a boundary condi- 
tion that depends on the outflow (instead of 5'|p = 0). Consistent with the 
comments of Lekien et al. 2004| , however, we found that the algorithm led to 



problematic numerical instability. The reconstructed velocity field varied widely 
from frame to frame in an essentially unpredictable way, leading us to use the 
velocity potential function method instead. 
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Figure 11: Regions of varying local inflow and outflow in an example velocity 
field as seen from above. An arrow represents the velocity of each tracked 
particle. The solid rectangle, cutting through vortex centers, comprises an open 
boundary with large local inflow and outflow. The dashed rectangle, passing 
near vortex separatrices, comprises an open boundary with small local inflow 
and outflow. 
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5 Conclusions 



The post-processing technique described in [JSjaUows the multi-scale reconstruc- 
tion of a velocity field using particle tracks. It requires neither interpolation nor 
smoothing, but automatically eliminates outliers and reduces spurious, appar- 
ent compressibility by projecting the measurements onto a set of incompressible, 
two-dimensional basis modes. By building a mesh from the Delaunay triangula- 
tion of the particle locations themselves, it is possible to calculate the necessary 
spatial gradients via finite-element techniques, without constructing an arbi- 
trary grid. A set of boundary modes 4>j approximates the velocity potential 
function, accounting for local inflows and outflows at the boundary. They are 
constructed numerically by solving Laplace's equation. The streamfunction, 
which accounts for the remaining flow, is expanded in a set of interior modes 
tpk- For these, the eigenmodes of the Laplacian form a natural basis. In ir- 
regular measurement regions, the eigenmodes may be calculated numerically. 
In regions of regular shape in which analytic eigenmodes are known, they may 
be used instead, reducing the computational expense. For our data, we have 
chosen Fourier modes, since the boundaries of our effective measurement region 
nearly match the rectangular field of view of our camera. 

The resulting reconstructed velocity fields give direct access to scale-local 
properties of the flow and can be decomposed without filtering. Spatial spectra 
can be constructed without windowing, boundary effects can be separated from 
the interior flow modes, and the technique holds promise for direct calculation 
of enstrophy and energy transfer rates. 

This technique can be applied to any set of two-dimensional flow measure- 
ments that can be expressed as a list of local velocity measurements at particular 
locations over time. In its current form, it is better suited for laboratory flows 
with high particle density than for oceanic flows with low density, since we 
have avoided supplementing the measurements with interpolation or modeling. 
Stratified flow experiments, rotating quasi-two-dimensional experiments, and 
soap film experiments are natural candidates for this type of post-processing, 
and we expect it to have wide applicability. 

This work was supported by the U.S. National Science Foundation under 
Grant No. DMR-0906245. 

References 

G. B. Arfken and H. J. Weber. Mathematical Methods for Physicists. Harcourt 
Academic Press, San Diego, fifth edition, 2001. 

C.-H. Bruneau, P. Fischer, and H. Kellay. The structures responsible for the 
inverse energy and the forward enstrophy cascades in two-dimensional tur- 
bulence. Europhys. Lett., 78(3):34002, 2007. URL http : //stac ks . iop . org/ 
10295-5075/78/34002, 



27 



E. BuUard and H. Gellman. Homogeneous dynamos and terrestrial magnetism. 
Phil. Trans. R. Soc. London A, 247(928) :213-278, November 1954. 

P. C. Chu, L. M. Ivanov, T. P. Korzhova, T. M. Margolina, and O. V. Mel- 
nichenko. Analysis of sparse and noisy ocean current data using flow decom- 
position. Part I: Theory. J. Atmos. Ocean. Tech., 20(4):478-491, 2003. 

H. J. H. Clercx, G. J. F. van Heijst, and M. L. Zoeteweij. Quasi-two-dimensional 
turbulence in shallow fluid layers: The role of bottom friction and fluid layer 
depth. Phys. Rev. E, 67(6):066303, Jun 2003. doi: 10.1103/PhysRevE.67. 
066303. 

R. H. Kraichnan. Incrtial ranges in two-dimensional turbulence. Phys. Fluids, 
10:1417-1423, July 1967. 

F. Lekien, C. Coulliette, R. Bank, and J. Marsden. Open-boundary modal 
analysis: Interpolation, extrapolation, and filtering. J. Geophys. Res., 109: 
C12004, 12 2004. URL http://dx.doi.org/10. 1029/2004JC002323, 

B. Liithi, A. Tsinober, and W. Kinzelbach. Lagrangian measurement of vorticity 
dynamics in turbulent flow. J. Fluid Mech., 528:87-118, 2005. 

P. Lynch. Partitioning the wind in a limited domain. Mon. Weather Rev., 117: 
1492, 1989. doi: 10.1175/1520-0493(1989)117. 

Manikandan Mathur, George Haller, Thomas Peacock, Jori E. Ruppert-Felsot, 
and Harry L. Swinney. Uncovering the Lagrangian skeleton of turbulence. 
Phys Rev. Lett, 98 (14): 144502 , 2007. doi: 10.1103/PhysRevLett.98.144502. 
URL |http : //link . aps . org/abstract/PRL/v98/el44502] 



N. Mordant, A. M. Crawford, and E. Bodenschatz. Experimental Lagrangian 
acceleration probability density function measurement. Physica D, 193:245- 
251, 2004. 

N. T. Ouellette. On the performance of predictive best-estimate particle- 
tracking algorithms, submitted, 2010. 

N. T. Ouellette, H. Xu, and E. Bodenschatz. A quantitative study of three- 
dimensional Lagrangian particle tracking algorithms. Exp. Fluids, 40:301-313, 
February 2006. doi: 10.1007/s00348-005-0068-7. 

N. T. Ouellette, H. Xu, and E. Bodenschatz. Measuring Lagrangian statistics in 
intense turbulence. In C. Tropea, A. L. Yarin, and J. F. Foss, editors. Springer 
Handbook of Experimental Fluid Mechanics. Springer- Verlag, Berlin, 2007. 

Nicholas T. Ouellette and J. P. GoUub. Curvature fields, topology, and the dy- 
namics of spatiotemporal chaos. Phys. Rev. Lett., 99(19):194502, 2007. doi: 
10.1103/PhysRevLett. 99.194502. URL http : //link . aps . org/abstract/| 
|PRL/v997el9 4502, 



28 



W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numeri- 
cal Recipes: The Art of Scientific Computing. Cambridge University Press, 
London, third edition, 2007. 

M. K. Rivera and R. E. Ecke. Pair dispersion and doubling time statistics in 
two-dimensional turbulence. Phys. Rev. Lett., 95(19):194503, Nov 2005. doi: 
10.1103/PhysRevLett. 95.194503. 

M. K. Rivera, W. B. Daniel, S. Y. Chen, and R. E. Ecke. Energy and enstrophy 
transfer in decaying two-dimensional turbulence. Phys. Rev. Lett, 90(10): 
104502, Mar 2003. doi: 10.1103/PhysRevLett.90.104502. 

L. Rossi, J. C. Vassilicos, and Y. Hardalupas. Electromagnetically controlled 
multi-scale flows. J. Fluid Mech., 558:207-242, July 2006. doi: 10.1017/ 
S0022112006009980. 

D. Rothstein, E. Henry, and J. P. GoUub. Persistent patterns in transient chaotic 
fluid mixing. Nature, 401:770-772, October 1999. 

T. H. Solomon and I. Mezic. Uniform resonant chaotic mixing in fluid flows. 
Nature, 425:376-380, September 2003. doi: 10.1038/nature01993. 

P. Tabeling, S. Burkhart, O. Cardoso, and H. Willaime. Experimental study of 
freely decaying two-dimensional turbulence. Phys. Rev. Lett., 67(27) :3772- 
3775, Dec 1991. doi: 10.1103/PhysRevLett.67.3772. 

F. Toschi and E. Bodenschatz. Lagrangian properties of particles in turbulence. 
Annu. Rev. Fluid Mech., 41:375-404, 2009. 

C. Tropea, F. Scarano, J. Westerweel, A. A. Cavone, J. F. Meyers, J. W. Lee, 
and R. Schodl. Particle-based techniques. In C. Tropea, J. Foss, and A. Yarin, 
editors. Springer Handbook of Experimental Fluid Mechanics, pages 287-362. 
Springer- Verlag, Berlin, 2007. 

Dominic Vella and L. Mahadevan. The "cheerios effect" . American Journal of 



Physics, 73(9):817-82 5, 2005. doi: 10.1119/1.1898523. URL |http://link. 
aip ■ org/link/?AJP/73/8177l| 



G. A. Voth, G. Haller, and J. P. GoUub. Experimental measurements of stretch- 
ing flelds in fluid mixing. Phys. Rev. Lett, 88(25):254501, June 2002. doi: 
10.1103/PhysRevLett. 88.254501. 



29 



